home *** CD-ROM | disk | FTP | other *** search
/ TPUG - Toronto PET Users Group / TPUG Users Group CD / TPUG Users Group CD.iso / AMIGA / AMICUS / AMICUS18.ADF / Progs / StarProbe / SPFIT.C < prev    next >
Text File  |  1989-01-27  |  5KB  |  227 lines
  1. /****************************************************************************
  2.  ***                                                                      ***
  3.  ***      S T A R   P R O B E        F I T T E R                          ***
  4.  ***                                                                      ***
  5.  ****************************************************************************/
  6. /*****
  7.  ***  This module evaluates how far the inward & outward integrations
  8.  ***  are from each other, and adjust the boundary values of (p,t,l,r)
  9.  ***  accordingly.
  10.  *****/
  11.  
  12. #include "stdio.h"
  13. #include "math.h"
  14. #include "limits.h"
  15. #include <dos.h>
  16. #include "spmac.h"
  17. #include "spref.h"
  18.  
  19. extern FILE *pf;
  20.  
  21. double delta     [4],
  22.        old_delta [4],
  23.        a         [4] [8];
  24.  
  25. /*****
  26.  ***  D i s c o n t i n u i t y
  27.  ***
  28.  ***  Returns (true) if any of the endpoints are more than <epsilon> apart.
  29.  *****/
  30.  
  31. int discontinuity()
  32. {
  33.  iam(70);
  34.  int i,count,val;
  35.  double epsilon;
  36.  block(in,"discontinuity",0.0);
  37.  
  38.  epsilon = 1.0e-2;
  39.  
  40.  count = 0;
  41.  for (i=0;i<4;i++){
  42.    old_delta[i] = delta[i];
  43.    delta[i] = results[i][0] - results[i][1];
  44.    if ((abs(delta[i]/results[i][0])) > epsilon)
  45.      count++;
  46.    }
  47.  
  48.  block(mid,"discontinuity count",(double)count);
  49.  val =  (count > 0) ? 1 : 0;
  50.  
  51.  block(out,"discontinuity",(double)val);
  52.  return(val);
  53. }
  54.  
  55. /*****
  56.  ***  X f e r   R e s u l t s
  57.  ***
  58.  ***  Transfers results array into 'a[]' 
  59.  *****/
  60.  
  61. void xfer_results()
  62. {
  63.  iam(71);
  64.  register int i;
  65.  double e1,e2,e3,e4;
  66.  block(in,"xfer results",0.0);
  67.  
  68.  e1 = corepress    * (adj_p - 1.0);
  69.  e2 = coretemp     * (adj_t - 1.0);
  70.  e3 = lumratio*SL  * (adj_l - 1.0);
  71.  e4 = radratio*SR  * (adj_r - 1.0);
  72.  
  73.  for (i=0;i<4;i++){
  74.    a[i][0] = (results[i][2] - results[i][0])/e1;
  75.    a[i][1] = (results[i][3] - results[i][0])/e2;
  76.    a[i][2] = (results[i][4] - results[i][1])/e3;
  77.    a[i][3] = (results[i][5] - results[i][1])/e4;
  78.    }
  79.  
  80.  block(out,"xfer results",0.0);
  81. }
  82.  
  83. /*****
  84.  ***  I n v e r t
  85.  ***
  86.  ***  Inverts the 4x4 matrix a[] into a'[]
  87.  *****/
  88.  
  89. void dump_a()
  90. {
  91. int i,j;
  92.  fprintf(pf,"\naaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa");
  93.  for (i=0;i<4;i++){
  94.    fprintf(pf,"\na");
  95.    for (j=0;j<8;j++)
  96.      fprintf(pf," %g",a[i][j]);
  97.    }
  98.  fprintf(pf,"\naaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa\n");
  99. }
  100.  
  101. void invert()
  102. {
  103.  iam(72);
  104.  register int i,j,im,ip,ist,n,n1;
  105.  double al;
  106.  block(in,"invert",0.0);
  107.  
  108. for (i=0;i<4;i++)
  109.    for (j=4;j<8;j++)
  110.      a[i][j] = (i == (j-4)) ? 1.0 : 0.0;
  111.  
  112.  if (trace_me) dump_a();
  113.  
  114.  n1 = 8; n = 4;
  115.  
  116.  for (ip=0;ip<n;ip++){
  117.    im = ip;
  118.    ist = ip + 1;
  119.    if (ip != n)
  120.      for (i=ist;i<n;i++)
  121.        if (abs(a[im][ip]) < abs(a[i][ip]))
  122.          im = i;
  123.    if (a[im][ip] == 0.0)
  124.      assert_error(proc_num,(double)im,(double)ip,(double)ist,0.0,0.0,0.0,0.0);
  125.    if (im != ip)
  126.      for (j=ip;j<n1;j++){
  127.        al=a[ip][j]; a[ip][j]=a[im][j]; a[im][j]=al;}
  128.    al = a[ip][ip];
  129.    a[ip][ip] = 1.0;
  130.    for (j=ist;j<n1;j++)
  131.      a[ip][j] /= al;
  132.    for (i=0;i<n;i++)
  133.      if (i!=ip){
  134.        al = a[i][ip];
  135.        for (j=ip;j<n1;j++)
  136.          a[i][j] -= al*a[ip][j];
  137.        }
  138.    }
  139.  if (trace_me) dump_a();
  140.  block(out,"invert",0.0);
  141. }
  142.  
  143.  
  144. /*****
  145.  ***  N e w   B o u n d a r y   C o n d i t i o n s
  146.  ***
  147.  ***  Ref: Sears & Brownlee (Aller) pp. 597-601
  148.  *****/
  149.  
  150. void new_boundary_conditions()
  151. {
  152.  iam(73);
  153.  register int i,j,k;
  154.  double new_val[4],dp,dt,dl,dr;
  155.  block(in,"new bdry cond",0.0);
  156.  
  157.  xfer_results();
  158.  
  159.  invert();
  160.  
  161.  for (i=0;i<4;i++){
  162.   new_val[i] = 0.0;
  163.   for (j=4,k=0;j<8;j++,k++)
  164.     new_val[i] += (a[i][j] * delta[k]);
  165.   }
  166.  
  167.  dp = new_val[0]/corepress;
  168.  dt = new_val[1]/coretemp;
  169.  dl = new_val[2]/(SL*lumratio);
  170.  dr = new_val[3]/(SR*radratio);
  171.  
  172.  block(mid,"new bdry cond +p",new_val[0]);
  173.  block(mid,"new bdry cond +t",new_val[1]);
  174.  block(mid,"new bdry cond +l",new_val[2]);
  175.  block(mid,"new bdry cond +r",new_val[3]);
  176.  
  177.  printf("... new bdry cond: %g %g %g %g\n",dp,dt,dl,dr);
  178.  fprintf(pf,"... new bdry cond: %g %g %g %g\n",dp,dt,dl,dr);
  179.  
  180.  corepress += new_val[0];
  181.  coretemp  += new_val[1];
  182.  lumratio  = ((SL*lumratio)+new_val[2])/SL;
  183.  radratio  = ((SR*radratio)+new_val[3])/SR;
  184.  
  185.  block(out,"new bdry cond",0.0);
  186. }
  187.  
  188.  
  189. /*****
  190.  ***  C l o s e r   F i t
  191.  ***
  192.  ***  Returns (true) if all (p,t,l,r) deltas were reduced from the previous
  193.  ***  base integrations (0 & 1).
  194.  *****/
  195.  
  196. int closer_fit(ixstop)
  197. int ixstop;
  198. {
  199.  iam(74);
  200.  int val;
  201.  double dP,dT,dL,dR;
  202.  block(in,"closer fit",0.0);
  203.  
  204.  dP = abs(pressgrad[ixstop]-pressgrad[ixstop+1])/
  205.       ((pressgrad[ixstop]+pressgrad[ixstop+1])*0.5);
  206.  dT = abs(tempgrad[ixstop]-tempgrad[ixstop+1])/
  207.       ((tempgrad[ixstop]+tempgrad[ixstop+1])*0.5);
  208.  dL = abs(lumgrad[ixstop]-lumgrad[ixstop+1])/
  209.       ((lumgrad[ixstop]+lumgrad[ixstop+1])*0.5);
  210.  dR = abs(radiusgrad[ixstop]-radiusgrad[ixstop+1])/
  211.       ((radiusgrad[ixstop]+radiusgrad[ixstop+1])*0.5);
  212.  
  213.  old_diff = cur_diff;
  214.  cur_diff = (dP+dT+dL+dR)*0.25;
  215.  
  216.  block(mid,"closer fit dP",dP);
  217.  block(mid,"closer fit dT",dT);
  218.  block(mid,"closer fit dL",dL);
  219.  block(mid,"closer fit dR",dR);
  220.  
  221.  val = (cur_diff < old_diff) ? 1 : 0;
  222.  
  223.  block(out,"closer fit",(double)val);
  224.  return(val);
  225. }
  226.  
  227.